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Abstract 

We present a numerical scheme for determining hyperboloidal initial 
data sets for the conformal field equations by using pseudo-spectral meth- 
ods. This problem is split into two parts. The first step is the determina- 
tion of a suitable conformal factor which transforms from an initial data 
set in physical space-time to a hyperboloidal hypersurface in the ambient 
conformal manifold. This is achieved by solving the Yamabe equation, a 
non-linear second order equation. The second step is a division by the 
conformal factor of certain fields which vanish on J ^ the zero set of the 
conformal factor. The challenge there is to numerically obtain a smooth 
quotient. Both parts are treated by pseudo-spectral methods. The non- 
linear equation is solved iteratively while the division problem is treated 
by transforming the problem to the coefficient space, solving it there by 
the QR-factorisation of a suitable matrix, and then transforming back. 
These hyperboloidal initial data can be used to generate general relativis- 
tic space-times by evolution with the conformal field equations. 



1 Introduction 

In this article we shall discuss the problem of numerically calculating "hyper- 
boloidal initial data" . These occur naturally in the numerical solution of the 
conformal Einstein equations, a promising approach towards the numerical evo- 
lution of general relativistic space-times ||7|, p| p^. 

Consider an asymptotically flat solution of the Einstein equation. Let S be 
a space-like hypersurface which extends to infinity in such a way that it reaches 
null infinity J . Assume that it remains space-like even in the limit then it 
will touch J transversely. Prototypes for such hypersurfaces are the space-like 
hyperboloids in Minkowski space, hence they are called hyperboloidal hyper- 
surfaces. These are rather different from the standard asymptotically euclidean 
hypersurfaces which end up at space-like infinity i". In contrast to the latter 
hypersurfaces the hyperboloidal hypersurfaces are not Cauchy hypersurfaces of 
the complete space-time for the standard Einstein evolution. However, they can 
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still be used for pre- or retrodiction depending on whether they end up on J^'^ 
or J-. 

The main reason for focusing on hyperboloidal hypersurfaces is the fact that 
outgoing radiation can be traced much more effectively on those hypersurfaces 
than it is possible on the asymptotically flat hypersurfaces ||T^. This is due 
to the fact that they approximate a null hypersurface at large distances along 
which the radiation propagates instantaneously. 

The initial data for the Einstein evolution on a hyperboloidal hypersurface 
S implied by the space-time geometry consists of the intrinsic metric and the 
extrinsic curvature. The requirement that S should reach out to null infinity 
imposes quite definite asymptotic fall-off conditions on these fields on S. After a 
suitable conformal compactification of the space-time the fall-off conditions can 
conveniently be captured in smoothness conditions for appropriately rescaled 
fields on the boundary of a three-dimensional manifold. In this "conformal" 
picture the evolution of initial data is most appropriately performed with the 
regular conformal field equations ^, |lO), yielding the "conformal method" for 
solving Einstein's equations. The main advantage of using this method is that 
the conformal field equations allow the evolution of initial data all the way 
up to time-like infinity , which is a regular point of the conformal manifold 
provided the data are sufficiently close to Minkowski data. But even if i+ does 
not exist one can follow the evolution up to the point where singularities form. 
This has been demonstrated in Thus, it is possible to evolve the complete 
future of an initial hypersurface on a single grid including the points of (and 
even beyond, since one can smoothly extend the initial data across J^). This 
allows physically meaningful quantities like the Bondi-news and -momentum, 
the radiation flux etc. which well defined only on to be determined there 
without any further approximation, essentially be reading off the appropriate 
functions from the grid. 

Of course, for the conformal field equations, initial data have to be deter- 
mined, too. They consist no longer only of the first and second fundamental 
form, satisfying the standard constraint equations. The conformal field equa- 
tions comprise more variables than the standard Einstein equations. In partic- 
ular, they include the Bianchi identity from which follow evolution equations 
for both the Weyl- and the Ricci curvature. Hence, an initial data set for the 
evolution with the conformal field equations is larger than for the standard evo- 
lution. Such an initial data set is called hyperboloidal if the initial surface is a 
hyperboloidal hypersurface. 

The conformal constraint equations, i.e., that part of the conformal field 
equations which is intrinsic to the initial hypersurface, need to be solved in or- 
der to obtain the initial data. At the moment there is no way to solve those 
constraints directly in more than one space dimension but there is a procedure 
due to Andersson, Chrusciel and Friedrich |^ which allows the construction of 
hyperboloidal initial data. Essentially, one solves one non-linear second order 
equation for one scalar function to determine a conformal factor which trans- 
forms from a solution of the constraints in physical space to unphysical space. 
Once this has been found one can determine the remaining hyperboloidal initial 
data algebraically from a part of the conformal constraint equations. 

In this procedure there are two complications. In the first place, the second 
order equation is singular at the boundary where the hyperboloidal surface 
intersects J^. Fortunately, this turns out to be more trouble for the analytical 
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than for the numerical treatment. And in the second place, when computing 
the initial data one has to divide by the conformal factor which vanishes on 

Although it can be shown analytically that this "division by zero" is well 
defined it does pose numerical problems. 

The purpose of this article is to discuss a numerical implementation for 
finding hyperboloidal initial data by pseudo-spectral methods and to suggest 
one possible way to overcome the division problem. In section ^ we describe 
the analytical background in more detail. The numerical solution of the second 
order equation is discussed in section |^ while the division problem is treated in 
section ^. A short discussion concludes the paper. 



2 Finding hyperboloidal initial data sets 

Let (M, g) be an asymptotically flat solution of the vacuum Einstein equation 
which has a smooth conformal structure at null infinity J^. Consider S, a 
hyperboloidal hypersurface intersecting in a smooth two-dimensional surface 
VS. The Lorentz metric g on M induces a riemannian metric h on E. Let 
be a conformal factor, which smoothly attaches null infinity to M thus defining 
the "conformal" space M = M UVM. Then it follows that E E U VS is 
a smooth submanifold with boundary in M and that > on E and = 0, 
dfl 7^ on VE. Furthermore, there exists a smooth riemannian metric ft, on E 
such that on E the relation 

h^n^h (2.1) 

holds. The embedding of E in M defines the second fundamental form fc on E. 
Being induced from a vacuum solution of the Einstein equation the pair [h, k) 
satisfies the physical constraint equations on E. 

When solving the constraint equations on asymptotically euclidean initial 
surfaces it is convenient to make the assumption of time-symmetry, namely 
that the extrinsic curvature of the initial surface should vanish. This condition 
is not compatible with the geometry of a hyperboloidal surface but one can 
obtain similar simplications in this case also by the assumption that the extrinsic 
curvature be pure trace, i.e., proportional to the metric, 

k = ^Kh. (2.2) 

Then the momentum constraint implies that if is a constant and therefore, so 
is the scalar curvature of h by virtue of the hamiltonian constraint. One can 
assume that (after rescaling h with a suitable constant) 
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i? = -6. (2.3) 



Note, that the assumption (2.2) which asserts that the tracefree part of the 
extrinsic curvature vanishes is conformally invariant so that the extrinsic cur- 
vature of E in the unphysical space is also pure trace, albeit not necessarily 
constant. 



In this paper we will impose the condition (2.2). This is yields simplest 
case for constructing hyperboloidal initial data sets. The analytic treatment 
of this problem has ben thoroughly discussed in B. But there are also several 
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other, less restrictive, treatments in the hterature. In Q the assumption ( |2.2| ) is 
dropped allowing for an extrinsic curvature which is almost general apart from 
the fact that the mean curvature is required to be constant. In also this 
requirement is dropped and in the existence of hyperboloidal initial data is 
discussed for situations with a non- vanishing cosmological constant. 

In order to construct hyperboloidal initial data one may proceed as follows 
[0 ^: let be a boundary defining function for E, i.e., > on S and 
w = 0, d^T^OonVE and let h he a smooth metric on E. Now we s eek a 



conformal factor so that the metric h = il^^h defined on S sat isfies (2.3). We 
write ri = (f)~'^uj for some smooth function ip on E. Then ( ^.3|) turns into the 
non-linear second order equation, also called the Yamabe equation 

4Lj^A<j) - 4wV°cjVa(/) - [cJ^R + 2wAw - BVcjVaw] (t> = 3(/>^ (2.4) 

Here, we have used the Laplace operator A = V^Va with respect to h, the co- 
variant derivative operator V and the scalar curvature R of h. The most obvious 
property of this equation is the fact that it degenerates on the boundary where 
u vanishes. This is not entirely surprising considering its origin: if the equation 
were regular on the boundary one would presumably be able to specify bound- 
ary data, thus introducing some freedom into the structure of a hyperboloidal 
hypersurface at its intersection with J^. But this would be in contradiction 
with the fact that is universal, being fixed entirely by t he smooth conformal 



structure. The conformal transfomation properties of (2.4) ensure that the con- 
formal factor defined from a solution (j) does not depend on the specific form 
of the boundary defining function uj and, furthermore, that it depends only on 
the conformal class of the metric h. 

Consider now the following fields defined on E 

= -1^"' (^VaVfef] - i/iafcAf^) , (2.5) 

Eab = (^Rab - ^habR - *a6^ • (2.6) 

These fields are the essential initial data necessary for the evolution with the 
conformal field equations, ^ab is the projection of the conformal Ricci tensor 
onto E, while Eab represents the rescaled electric part of the Weyl tensor. As 
they stand these expressions are valid only on E, being formally singular on 
the boundary where f2 vanishes and one needs to worry whether there exists a 
smooth extension to VE. This question and the more immediate question of 



existence, uniqueness and regularity of solutions of (2.4) have been answered in 



complete detail in |g| where the following theorem has been proved. 

Theorem 2.1 Suppose (E,/i) is a three-dimensional, orientable, compact, smooth 
Riemannian manifold with boundary VE. Then there exists a unique solution 



of (2.4) and the following conditions are equivalent: 



1. The function (p as well as the tensor fields (2.5) and ( |2.6| ) determined on 
E from h and — (j)~'^uj extend smoothly to all of E. 

2. The conformal Weyl tensor Cab = ^Eab goes to zero on VE. 

3. The conformal class of h is such that the extrinsic curvature of VE with 
respect to its embedding in (E, h) is pure trace. 
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Condition 3. is a weak restriction of the conformal class of the metric ft, on S, 
since it is only effective on the boundary. Interestingly, the theorem only requires 
S to be orientable and does not restrict the topology of S any further. We exploit 
this fact by assuming the existence of an isometric and hypersurface orthogonal 
action of C/(l) without fixed points on E which we take to have topology x 

X /. Thus, we may ignore the dependence on one coordinate, reducing the 
problem to one on the two-dimensional surface x /. We take coordinates 
on this surface as {u,v) G [0, 27r] x [—1,1], with 27r-periodicity implied on the 
u-dependence. The boundary is given hy v — ±1. In the case treated here in 
detail, we choose J to coincide with the boundary. We refer to Q for a more 
detailed discussion of space-times with these properties. In more complex cases 
one could think of choosing boundary functions u which vanish not only on the 
boundary but also somewhere in the interior. This would define and evolve to 
physically interesting space-times with more complex geometries We will 
show one possibility below. 

We choose the boundary function to be 

'^{u,v) = i (1 . 



On the boundary the equation (2.4) reduces to 

(V"wVaw)0 = 0^ (2.7) 

Since (j) has to be non-zero on the boundary this implies (j) = ^VwVaW, so 
that the boundary values of the solution are completely fixed by the equation. 
Due to our assumptions the metric on E has the form 

h = huu du^ + 2huv du dv + h^y dv^ + h^^ dw'^ , 

where the metric functions do not depend on the coordinate w. Since we need to 
specify only the conformal class of a metric on E we may assume that h^^j^^ — 1, 
leaving the other functions arbitrary except for the condition 3. of Theorem 



2.1 . With these assumptions, the induced metric on the boundary is 

P ~ huu du^ + dviP' , 

and the extrinsic curvature of the boundary with respect to the metric h is 
proportional to 

A — f/?™/) +h'"''h +2h ft™ ) dy^ 

Condition 3. requires that A be proportional to the induced metric p which 
implies, that A itself has to vanish. One possibility to satisfy this condition is to 
require that h^y — and huu.v = on the boundary. It is worthwhile to point 
out again that the purpose of condition 3 is to ensure the smooth extensibility 



of the solution and the tensor fields mentioned in theorem 2.1. It is possible to 



find solutions of (2.4) for free data which do not satisfy condition 3. 



3 Numerical solution of the Yamabe equation 

We have implemented a numerical scheme for the solution of ( |2.4| ) based on 
pseudo-spectral methods. 
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In this section we want to describe a numerical scheme based on pseudo- 
spectral methods for solving the Yamabe equation. There are various reasons 
for considering these methods. They are known for their high accuracy, at 
least in situations where the solution is smooth. Then the numerical error 
decreases exponentially with the number of degrees of freedom (i.e., the number 
of collocation points or the number of basis functions used to approximate the 
solution). This is much faster than the error decay in any finite difference 
method which is 0{N~'^) with q usually less than 4. This property makes 
pseudo-spectral methods ideally suited for elliptic problems. Pseudospectral 
methods have been employed successfully in various areas of physics and applied 
mathematics. In particular, we mention the work of S. Bonazzola and his co- 
workers e.g., Q on various applications in relativistic astrophysics. 

Let us briefly describe the basic idea behind the use of spectral and pseudo- 
spectral methods. For more information on these methods we refer to the stan- 
dard textbooks ^, When solving a partial differential equation (time 
independent for our immediate purposes) in some domain S 



where C is some non-linear differential operator one seeks a solution in the form 
of an expansion in some suitable basis functions (assumed to be a complete set 
on the region of interest) 



The most popular functions in use are the trigonometric polynomials e""^ and 
the Chebyshev polynomials Tm{x). The choice depends on the topology of 
the domain and the boundary conditions. Inserting this expansion into the 
PDE yields a system of equations for the expansion coefficients. There are 
various ways to set up these equations. Spectral methods such as the Galerkin 
method reexpand ^iJ2n=i fk4>k{x)) in terms of the basis functions. This is 
only realistic in a few cases, mostly if C is constant and linear. In most other 
cases, the determination of the expansion coefficients of £f in terms of those 
of / is either impossible or computationally too expensive. Then one can fall 
back on the pseudo-spectral or collocation method: one introduces suitable 
collocation points xi, . . . ,xn and then the approximate solution is forced to 
satisfy the equation at the inner points and the boundary conditions. This 
yields N equations for the N expansion coefficients. 

The existence of the collocation points allows dual representations of the 
function /. Besides the "physical" representation based on the function values 
f{xi) at the collocation points there is the "spectral" representation based on the 
expansion into basis functions. The idea of the collocation method is to switch 
freely between those two representations using whichever is best to evaluate the 
various terms in the operator. This is made possible (at least for the Fourier 
and Chebyshev polynomials) by fast Fourier transformation (FFT) techniques. 

The representation in coefficient space is ideally suited for efficiently and 
accurately evaluating derivatives. The coefficients of the derivative of a func- 
tion are easily determined from the coefficients of the function either by simple 
multiplications or else by three term recurrence relations. Nonlinear and/or 



Cf^O 



(3.1) 



N 




(3.2) 
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nonconstant terms in the operator are best computed using the physical repre- 
sentation. 

In general, the matrices which represent the spectral approximations of even 
the simplest linear operators are full and difficult to invert directly. Their effi- 
cient inversion can be achieved by choosing a suitable preconditioning operator 
(see below). 

The method we employ for solving the Yamabe equation is an iterative 
method using a defect correction scheme. It is based on the following observa- 



tions. We write equation (2^) formally as 

C(j)^(j)^. (3.3) 

Here, >C is a linear operator made up from derivative and multiplication opera- 
tors. We construct a Richardson iteration procedure by writing (f>"^^ = 0" -I- (50. 



Inserting this into (3.3) and ignoring terms of higher order than the first in the 



correction term S(f) yields 

£50-5(0")''(50=-(/:0"-(0")^). (3.4) 

Thus, the general procedure for solving ( ^.4| ) is the following. Suppose we have 
some suitable approximation and compute the residual r" — C(f>^ — ((/)") . 
Then solve the linear equation ^£ — 5 (0")*^ S(j) = — r" to obtain the correction 

and an updated guess 0"+^ = -|- Scj). With that repeat the procedure until 
the accuracy is satisfactory. We observe that the linear operator acting on (50 
is also updated at each step but only by diagonal terms. 

As pointed out above, for pseudo-spectral methods, the matrix representa- 
tion of the operator £ is generally a full N x N matrix Cps, so that inversion 
is prohibitive both in terms of time and storage requirements for high N and 
especially in higher dimensions. However, there exists a way to circumvent this 
problem which is due to S. Orszag |2^: instead of inverting the pseudo-spectral 
representation of £, we substitute a finite-difference approximation CpD of £ 



into the left-hand side of equation (3^) and use this for the iteration procedure. 
In general, FD-approximations have sparse matrix representations, so that the 
iteration equation can be solved efRciently. Note, that this substitution is only 
made on the left-hand side of the equation while on the right-hand side the full 
pseudo-spectral approximation is retained. As pointed out in [ p2[ this method 
allows the efficient solution of general problems with operation costs and storage 
not much larger than those of the simplest finite difference approximations to 
the problem with the same number of degrees of freedom. 

In a sense this is similar to an inexact Newton method where the exact 
Jacobian is replaced by some approximation. Under such circumstances one 
cannot expect to have the full quadratic convergence of the Newton method, 
for which successive errors satisfy e„+i < Ke^ for some K > Q. Instead, under 
most circumstances one can hope for linear convergence, i.e., e„+i < Kcn, so 
that e„ < q;" for some positive a < 1 |Q . 

The solution procedure outlined in the previous section is implemented as 
follows. The topology suggests that we expand the fields into a Fourier series 
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in the periodic u-direction and Chebyshev polynomials in the w-direction as in 



M 



Tm{v). 



(3.5) 



We introduce the collocation points {ui, Vk), where Ui = ^i, i = 0, . . . ,N — 1 
and Vk = cos(^), fc = 0, ... ,M. Then we can switch between the physical 
and the spectral representation by using fast transformation techniques. The 
free data huu, h^v, h^v are specified on the collocation points subject only to 
the conditions huv = and huu,v = on the boundary. From the metric 
functions the connection (i.e., the Christoffel symbols) and the scalar curvature 
are determined in order to obtain the differential operator C. As indicated 
above, we use both its spectral approximation as well as a finite difference 
approximation. In the present case, we take an approximation which is 2"''- 
order accurate. In deriving the expression for this approximation one has to 
take into account that the collocation points Vk are not uniformly distributed. 

The equation (3.4) is then imposed at all the interior collocation points ex- 
cept at the boundary, where the values for the solution determined from (2.7) 
are inserted. This yields a matrix of size {N{M — 1))^, which is sparse. The 
computation of the residual, i.e., the right hand side of ( |3.4|) at each iteration 
step is done with the full spectral accuracy. The linear equation for the cor- 
rection is solved iteratively by methods taken from the sparse matrix package 
LASPACK [ p5| 

In Fig. hi is shown the convergence behaviour for a a typical run with 32 




100 200 300 

No. of Iterotions 



Figure 1: Size of residual and of update in the oviter loop versus the number of 
iterations. 

degrees of freedom (i.e., base functions) in each coordinate direction. The upper 
line is the logarithm (base 10) of the maximum size of the residual, while the 
lower line shows the logarithm of the maximum size of the update at each 
iteration step of the outer loop. Consistently, the residual is about three orders 
of magnitude above the update. The convergence is exponential according to 
the formula 



A cx a 



N 



(3.6) 
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where we find that in this case a « 0.96. The value of a depends on the number 
of degrees of freedom M = 32 • 33 and the free data. As mentioned above, the 
convergence behaviour is not the one expected for a true Newton method. This 
is not too surprising since the hnear operator we use for obtaining the update 
at each iteration step is not the Jacobian of the non-hnear function. When the 
correction hits the level of machine accuracy no further reduction in the residual 
is possible and the convergence levels off, the residual remaining at a level of 
about 10~^^. The exact numbers depend on the free data. 

In Fig. ^ are shown two solutions of the Yamabe equation for different kinds 
of free data. For the left diagram (case 1) we have chosen the following data 





Figure 2: The conformal factor obtained from the same expres- 
sions for the metric functions but different boundary functions. 
The region > corresponds to the physical space-time. 

huu = 2 + cj^ (v^ - sin^ u)^^ , 
huv^2uj (v^-sin^u) e~"'=°''", 

L n ^ — "^v cos u 

The right diagram (case 2) is meant to be a demonstration of the possibility to 
specify boundary functions which also vanish inside the computational grid. It 
was obtained by choosing 




as the boundary function and keeping the same expressions for the metric func- 
tions. It is apparent from the contour plots that the different boundary functions 
correspond to different topologies of the interior regions. For case 1 there are 
two J7's at the grid boundaries v = ±1, so that the physical space-time has 
the topology x /, as intended. But in case 2, there are additional zeros of 
w which introduce another J' inside the grid boundaries which has a circular 
topology. The resulting conformal factor defines the physical region inside that 
circular J^, giving it the topology of a disc. This two-dimensional picture is then 
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swept around by the J7(l)-action to yield a three-dimensional space-time with 
the topology of a full torus whose boundary is J^. This physical space-time is 
embedded in an "unphysical" region which is again bounded hy a on each 
side. 

We keep the i/'s at the grid boundaries because they prevent us from having 
to specify boundary condtions. Choosing the boundary function appropriately 
can lead to rather complicated space-times: by simply changing the sign of the 
boundary function in the example above we can switch the interior (physical) 
and exterior (unphysical) regions. In this way we obtain a space-time with 
toroidal J''s and another "asymptotically flat" end, which one could loosely 
interpret as a black hole region. However, one should be careful with inter- 
pretations like this because we still have the assumption of the existence of a 
hypersurface orthogonal Killing vector. 



4 The division problem 

We now turn to the problem of calculating the remaining initial data from a 
solution of the Yamabe equation. As indicated in section || this involves a 



division of the tensor components (|2.5|) and (2.6) by the conformal factor. From 



theorem 2.1 one knows that these components vanish on J' and that the quotient 
is smooth across J^. Let us first discuss a one-dimensional example. 

Consider two real valued functions on the interval [—1, 1] with /(O) — g{Q) — 
and f{x) ^ 0, g{x) ^ elsewhere. Our task is to compute the quotient f /g 
on the interval. The only problem is at the origin a; = 0, where the quotient 
is not defined. Analytically, one can use I'Hopital's rule to obtain the limit 
lmix^o{f / g){x). Numerically, however, the problem is more subtle, at least if 
one is interested in obtaining an answer as accurately and smoothly as possible. 

Let us be more specific in our assumptions on / and g. We take them both 
to be at least and to vanish at x = 0, but with g'{0) ^ 0. Then we have 
/(x) = xf(x) and g{x) — xg{x) for C^-functions / and g with g(0) ^ 0. Then 

lim(//5)(a:) = 

but obviously this limit cannot be calculated numerically in a direct way. The 
limit procedure has to be realized somehow. A straightforward method would 
be to approximate 

5(0) g(e)-5(0) g{e) 

for small enough e. Unfortunately, this is only a first order approximation which, 
of course, could be improved by using more accurate finite difference formulae 
for approximating the derivatives at a; = 0. Still, we do not get the accuracy of 
a spectral method. 

We propose here a method which is more in line with the idea of spectral 
methods. Roughly speaking, we transform the problem to coefficient space, solve 
it there and then transform back to physical space. Define Tm =^ span(T„j, < 
m < M), the space of polynomials on [—1,1] with degree at most M. This 
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space is spanned by the first AI + 1 Chebyshev polynomials and it carries a 
scalar product defined by 

/■I dx 

{Tni\Tn) = / Tn{x)Tm{x) 



The Chebyshev polynomials are orthogonal with respect to this scalar product, 
but not normalized. For each g G Tm, multiplication by g defines a linear map 
g : Tm T2m- This map and its matrix representation with respect to the 
basis polynomials follow from the Clebsch-Gordan like formula 

Thus, e.g., the matrix representation of (multiplication by) Tq is the (2 A/ + 
1) X [M 4- l)-matrix which is the identity in the upper half and zero otherwise, 
while Ti and T2 have the following (2M + 1) x (M + l)-matrix representations, 
respectively 



/ 1 \ 

2 1 
1 1 
1 1 



V 



/ 1 \ 

1 1 

2 1 
1 1 

V ■ • ■ • / 



(4.2) 



The higher degree polynomials T„ have similar representations. A general poly- 
nomial g G Tm is a linear combination in the Chebyshev polynomials and hence 
its matrix representation G is obtained by adding up these basic representations 
appropriately. Having the representation with respect to the basis polynomials 
it is easy to obtain the representation with respect to the normalised Chebyshev 
bases in Tm and T2m- 

Suppose now that / is in T2m- Obviously, the image of Tm under multipli- 
cation by g is an {M + l)-dimensional subspace of T2M, hence not all / G T2M 
are also in that image. Our task is to invert the map g : Tm ^[Xa/] on its 
image. Thought of in terms of linear algebra, this requires to solve (2M -I- 1) 
linear equations, of which only (M -I- 1) are linearly independent for (M -I- 1) 
unknowns. The matrix G is the coefficient matrix of that system of equations. 

We can solve this problem by finding the reduced QR-factorisation of G, 
i.e., we seek matrices Q and R so that that G = QR, where R is an upper 
triangular (M -I- 1) x (Af -I- l)-matrix and Q is a (2Af + 1) x (Af -I- l)-matrix 
whose columns are mutually orthonormal. Thus, the columns of Q form an 
orthonormal basis for the image of G. It should be noted that the scalar product 
involved here is the one defined between the Chebyshev polynomials if one uses 
the normalised polynomials when representing matrices, which will be assumed 
in the sequel. For a more thorough discussion of the QR-factorisation from 
various perspectives we refer to standard textbooks on numerical linear algebra 
like 0, 0, |§. 

Suppose now that G — QR has been factored. We note that Q*Q = 1m, 
the identity in T/v/ , while QQ* is the orthogonal projector onto the image of 
G. Given any / g 72a/, the QR-factorisation enables us to "solve" the overde- 
termined system Gh = / in the following way. We have QQ*/ = Gh for some 
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h G Tm- Hence we get Q*/ = R/i and, finally, inverting R, we get h. It is a 
well known property of the QR-factorisation that it allows the solution of least 
squares problems, i.e., given the overdetermined equation Gh = /, the "solu- 
tion" h = R~^Q*/ has the property that it minimizes \\Gh — Thus, e.g., 
if f is in the image of G then h is the (unique, if R is invertible) vector in Tm 
for which the equation holds. But if / is not in the image, then h will be that 
vector in Tm whose image is closest to /, i.e., for which / — Gh is orthogonal to 
the image. Therefore, the QR-factorisation serves to compute the "generalized 
inverse" or Moore-Penrose inverse of a matrix p3| . 

Now we can find the solution of the division problem as follows. Given a 
smooth function / on the interval [—1,1] we compute its expansion into Cheby- 
shev polynomials up to degree M. Then we can consider / to be an element 
of Tm- We may also regard it as being in Tzm by taking the coefficients of the 
polynomials with degree higher than M to be zero. We also expand the divisor 
g, another smooth function on the above interval and from its expansion coef- 
ficients we construct the matrix G. Then we compute the QR-factorisation of 
G by standard methods (we take Householder reflections). The solution {f /g) 
is obtained by computing y — Q*/, solving Rj/ = z by back substitution and, 
finally, by transforming back from the expansion coefficients contained in the 
vector z to the function {f /g). 

Let us illustrate the above mentioned behaviour by an example. We take 
g{x) — X and f{x) = sin(10a:) -I- 2cos(5a;) — 2. Since /(O) = 0, / is divisible 
by g with {f/g)(0) = 5. In Fig. ^ is shown the exact quotient (solid line) 




Figure 3: The quotient (f/g) (see text) 

and three approximations for h = (f/g) with M = 8,16,32 obtained by the 
procedure given above. Obviously, the lowest approximation with M = 8 is 
not usable. The reason for this is that there is too much structure present in 
/ to be resolved with only 8 polynomials. Following the discussion in |2^ and 
[ p^ , we find that in this case we need about 10 polynomials to have enough 
resolution power, and indeed, the approximation with 16 polynomials is almost 
indistinguishable graphically from the exact function. With 32 polynomials the 
residual || / — /i • g||oo is on the level of the machine precision, see Table ^. To 
illustrate the behaviour when / is not in the image of g, we take f{x) = 1 and 
g{x) = X. The result is shown in Fig. ^ for M = 32. The thin line is the 
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M 


\\f -h-g\\oo 


8 


1.64 


16 


2.19(-3) 


32 


3.12(-14) 



Table 1: Maximum residual for approximations with different 
number M of polynomials. 




Figure 4: The quotient 1/x (see text) 

exact function 1/x while the thick line is the computed approximation h. The 
markers indicate the values at the collocation points Xk = — cos(fc7r/M), where 
the approximation is very good. However, inbctwccn the collocation points 
the agreement is bad because of the high-frequency oscillation. The maximal 
residual is in this case ||1 — x-/i||oo = 1-0, which is the value of / at the zero 
of g. But the maximum value of h is much higher, ||/i||cx3 = 10.52, one order of 
magnitude above the value of / at the zero of g. 

In this case the zero was in the center of the interval. If, instead, g vanishes at 
the boundary, then the singular behaviour is much more prominent. This is due 
to the clustering of the collocation points towards the end of the interval which 
normally is a benefit, because it allows a much more accurate approximation 
of functions there. In our case it means a good approximation of the singular 
behaviour. If we take f{x) = 1 but g{x) = 1 — x then we obtain a quotient with 
a maximum value of 373.9, more than two orders of magnitude above the value 
of / at the zero of g. 

If we compute (//<?) numerically from a function / which is known analyti- 
cally to have common zeros with g, but has been obtained by numerical means 
then / will in general not vanish exactly at the zeros of g. At those points its 
value e will depend on the accuracy of the algorithm used to compute /. Thus, 
we will have f{x) = g{x)h{x) + e (assuming for the moment that g has only one 
zero). Then, using the procedure described above, we would obtain h together 
with some additional "singular" or "high-frequency" part, which contaminates 
the result because its maximum value can be several orders of magnitude be- 
yond the value of e. The reason for this is that / is projected orthogonally onto 
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the image of g. But we are not so much interested in the orthogonal projection 
but rather we seek a projection which annihilates that high-frequency part. 

This can be achieved by a simple alteration of the projector P = QQ* in 
the following way. We seek a projector P which has the same image as P but 
which annihilates a vector, say y. We make the ansatz 

Pz = 'Pz - {a\z)Py, (4.3) 

where a is some vector with {a\y) — 1. For this to be a projector we need in 
addition that {a\Pz) = for all z. Thus, a is orthogonal to the image and the 
simplest choice for it is 

« = /^. (4.4) 

\\y - Pylr 

Inserting the expression for P = QQ*, one readily finds 



Pz = Q 



^ {y\y) - {Q*y\Q*y) ^ \ 



Now it is clear what to do numerically. First the vector y is determined as the 
expansion coefficients of 1 with respect to the normalised Chebyshev polyno- 
mials and, once the QR-factors of G have been found, the denominator of the 
factor above is computed. Then each dividend is Chebyshev expanded to get 
the vector z of its expansion coefficients, the term in brackets above is computed 
and, finally, the solution is obtained by inverting R. If this procedure is applied 
to find 1/x as before, then one obtains exactly zero. And in case one tries to 
find {f{x) + e)/g{x) the result is the same no matter what value of e. 

So far, this procedure works only for the case where g vanishes at a single 
point. If there are more zeroes of g present then the procedure has to be 
modified. This modification is straightforward. For two vectors yi, y2 to be 
annihilated, we have the altered projection 

Pz^Pz+-L {(yi\z)(y^\yi) ~ {yi\z)(yi\yi)) Py, 



2 



+ 7^ {{yi\^){yUyi) - (yil^Kyilyi)) Py2, (4.5) 

where V12 = {yi\yi) {y2\y2) - {y2\yif and where y^ = y -Py \s the part of 
y, orthogonal to the image of P. In principle, this can be generalized to even 
more y's but the formulae become more and more complicated. The meaning 
of the vectors yi is the following. Suppose g has n zeros in the interval. Then 
it is described by a polynomial of degree n. Each function in the image of G 
is then necessarily a polynomial of degree at least n. Therefore, the projection 
onto the image of G has to annihilate all the lower degree polynomials. This 
implies that one can take the standard basis vectors (1,0,...), (0, 1, ... ), ... 
for the vectors y. 

Let us now describe how to implement this method of division. We assume 
that a solution (/> of the Yamabe equation has been obtained. Together with 
the boundary defining function oj it defines the conformal factor f2 — ujcj)^'^. 
Since > 0, both O and to vanish at exactly the same points. From cf) and the 



geometry of S one determines the components of the tensor fields in (2.5) and 
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by differentiation and algebraic manipulations. Let ip be any one of those 
components. Analytically, one knows that V' shares the same zeroes as if the 
free data have been specified appropriately, so that, ideally, there is no problem 
when dividing ip by fl. However, the values of ip will never be exactly zero at 
the zeroes of $7 and we have the situation described above. 

We now assume for simplicity that a; is a function only of v, which is the 
case in the first example in section ^. This assumption is easily removed. To 
divide hy fl we divide by uj and then multiply with cf)'^. The function tp is 
represented by a two-dimensional array. To divide this array by lo we divide 
it row by row, each row consisting of the values Vi = ^(ui) at the points 
of constant u = Ui. The boundary function lu is Chebyshev-transformed and 
from its spectral representation we construct the matrix G and its reduced QR- 
factorisation. In the example with uj{v) = — w^) there are only two expansion 
coefficients because uj is an even quadratic polynomial. Then we compute the 
projection ( |4.5| ) of the Chebyshev-transformed row ipi onto the image of G. The 
two vectors yi and y2 have components Sko and 6ki respectively, thus being the 
spectral representation of the two lowest degree polynomials. The projection 
corresponds to the subtraction from 'ijji of an affine function av + b which agrees 
with ipi at the zeroes v — ±1 of uj. Thus, the projected function vanishes on the 
boundary and we can divide it by inverting the QR- factors of G. Performing 
these operations on all rows of tp finally yields the result ip/uj. 

As an example we show in Fig. ^ a component of the rescaled Weyl tensor 
which involves two divisions by lu, see eqns. (2.5) and (2.6). The maximal 
residual is in this case ||'0 — • a;||oo ~ 4(— 16) which is on the level of machine 
accuracy. This function is obtained from the first solution of the Yamabe equa- 




Figure 5: A component of the rescaled Weyl tensor 



tion (2.4) presented in the previous section. There are no interior j/'s and so 
the boundary of space-time coincides with the boundary of the grid. 



5 Conclusion 

We have presented in this paper a way for obtaining hyperboloidal initial data 
sets for the conformal Einstein equations by using pseudo-spectral methods. It 
has been demonstrated that these can be efhcient and powerful tools also in 
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numerical relativity. We have presented some results obtained under certain 
simplifying assumptions. In particular, these are the dimensional reduction 
by introducing a hypersurface orthogonal symmetry and the use of periodic 
boundary conditions. 

Although the results presented here are encouraging there are some points 
to be stressed. One desirable thing for the evolution of these initial data is 
that they be extended beyond J' in a, way which is as smooth as possible. This 
means that one has one or more J^'s in the interior of the grid just as in the 
second example presented in section ^ and then the boundary of the grid is 
outside the physical space-time. If the Yamabe equation is not singular at the 
grid boundary this means that we have to give some boundary values. These 
in turn influence the solution up to the inner J''s while inside these jT's the 
solution is determined from the equation alone. In general, these two parts of 
the solution will not match smoothly across the inner boundaries, there will be 
some higher derivative which jumps. This is due to the fact that the third one- 
sided normal derivative of the solution on J7 is characterised by the mass aspect 
of that part of the region on which the derivative is taken [|2j . And since there 
is no obvious relation between the values of the solution on the grid boundaries 
and the mass aspect on J' there is no guarantee that the third derivatives of 
the solution taken on either side of agree on J'. This implies that the initial 
data, and most notably the Weyl curvature, are not as smooth across as they 
should be. 

There might be a possibility, though, to overcome this problem (due to S. 
Brandt [Q). If we keep the grid boundaries singular then we do not need to 
specify boundary conditions there. Then it might be possible that the specifica- 
tion of free data which are analytic enforces the equality of the mass aspects of 
the different regions on their common J7. This possibility needs to be studied 
in detail. 

We have shown here some results in two dimensions. The generalisation 
of the Yamabe solver to full three dimensions should be straightforward. The 
situation is not so clear for the divison part. The fact that the two-dimensional 
division is performed by stacking together one-dimensional divisions might be 
impractical in three dimensions. 

Another limitation of the division procedure is the fact that it has to be 
changed whenever the boundary defining function acquires more zeros. This 
is more a matter of practicality than of principle. A final remark concerns 
the accuracy of the division method. Although each individual division can be 
made quite accurately with this method this does not imply that the quotients 
are similarly accurate, in the sense that the conformal constraint equations are 
satisfied to any accuracy comparable to that of the divisions. The reason for this 
is that the division process only "cures the symptoms" , namely the mild non- 
vanishing of the dividends on the boundary. It does not remove the reason for 
this phenomenon. That would probably be related to the fact that the Yamabe 
equation does not only enforce the behaviour of the function on the boundary 
but also of its first derivative. This has not been used so far in the solution 
process of the equation. Therefore, we cannot expect that the derivative of the 
numerical solution has the values it should have on the boundary. And so the 
constraints are only satisfied to the degree with which these values are attained. 
This problem is currently being studied. 

Thus, to summarize, we feel that the Yamabe solver can be made a rather 
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efficient and accurate tool, while the divisor may have its limitations. These 
come mostly from the fact that it cannot be easily applied for more general 
boundary functions. There is a completely different approach to the division 
problem developed by Hiibner p^ , where one solves a linear elliptic equation 
for the quotient ip/^^ which is singular at a; = 0. 

The results and methods presented in this paper demonstrate the feasibility 
of numerically determining initial data sets for the conformal field equations 
along the lines of II, |. They enable the numerical evolution of general rela- 
tivistic space-times using the "conformal method" . 
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